###
# R script for replicating empirical analysis in "Q-Monetary Transmission", by Priit Jeenas (priit.jeenas@upf.edu) and Ricardo Lagos (ricardo.lagos@nyu.edu). May 30, 2023.
###
library(plm)
library(plyr)
library(latex2exp)
setEPS()

# Define lag function that behaves as in earlier R versions
rlag <- function(x, k=1, shift="row", ...) {
  lag(x, k, shift, ...)
}
options(width=120)
# Pre-set parameters for figures
linew_main = 2.2; linew_ci = 1.6; col_change_ind = TRUE
font_ax_val = 1.5; font_leg_val = 1.4; font_lab_val = 1.5; linew_main = 3.0

# Reading data
data <- read.csv("JL_QMT_data.csv")

# Construct indicators for SIC-codes
if (TRUE){
    data$sic_str <- as.character(data$sic)
    # Add zero to short SICs
    data$sic_str[nchar(data$sic_str) < 4] <- paste("0", data$sic_str[nchar(data$sic_str) < 4], sep="")
    # Create broader SICs
    data$sic1 <- as.integer(substr(data$sic_str,1,1))
    data$sic2 <- as.integer(substr(data$sic_str,1,2))
    data$sic3 <- as.integer(substr(data$sic_str,1,3))

    data$sic <- as.factor(data$sic)
    data$sic1 <- as.factor(data$sic1)
    data$sic2 <- as.factor(data$sic2)
    data$sic3 <- as.factor(data$sic3)

    # Paste together the sic and datacqtr
    data$sic_datacqtr <- paste(data$sic_str, as.character(data$datacqtr), sep=",")
    data$sic_datacqtr <- as.factor(data$sic_datacqtr)
    data$sic1_datacqtr <- paste(as.character(data$sic1), as.character(data$datacqtr), sep=",")
    data$sic1_datacqtr <- as.factor(data$sic1_datacqtr)
    data$sic2_datacqtr <- paste(as.character(data$sic2), as.character(data$datacqtr), sep=",")
    data$sic2_datacqtr <- as.factor(data$sic2_datacqtr)
    data$sic3_datacqtr <- paste(as.character(data$sic3), as.character(data$datacqtr), sep=",")
    data$sic3_datacqtr <- as.factor(data$sic3_datacqtr)
}

###
# Construct new variables
###
# Recompute real assets based on GVA deflator
data$atq_real   <- 100 * data$atq / data$GVADef
# Drop negative capxq and sppeq observations
data[data$capxq<0 | is.na(data$capxq),"capxq"] <- NA
data[data$sppeq<0 | is.na(data$sppeq),"sppeq"] <- NA
# Same for capxq and saleq
data$capxq_real <- 100 * data$capxq / data$GVADef
data$sppeq_real <- 100 * data$sppeq / data$GVADef
# Leverage (right away in percentages)
data$lev_rat  <- 100*(data$dlcq + data$dlttq)/data$atq
# Cash-to-assets (right away in percentages)
data$cheat_rat <- 100*data$cheq/data$atq
# Log capital
# Real (log) sales
data$saleq_real <- 100 * data$saleq / data$GVADef
data$lsaleq_real <- log(data$saleq_real)
# Extra balance sheet variables
data$olq    <- data$ltq - data$dlcq - data$dlttq        # "Other liabilities" (other than debt)
# Take logs of balance sheet variables, and multiply right away with 100 (for balance sheet responses)
data$lppentq <- 100*log(data$ppentq)
data$latq   <- 100*log(data$atq)
data$lltq   <- 100*log(data$ltq)
data$ldttq  <- 100*log(data$dttq)
# Take other relevant ratios wrt assets (for balance sheet responses)
data$ltat_rat     <- 100* data$ltq / data$atq
data$olqat_rat    <- 100* data$olq / data$atq
data$ppentqat_rat <- 100* data$ppentq / data$atq

###
# PANEL CREATION
###
# Create panel dataframe only to compute certain lagged variables
p_data <- pdata.frame(data, index=c("gvkey", "cqtr_num"))

# Construct measures of equity issuance
if (TRUE){
    # First, compute gross issuances and repurchases
    p_data$eqissq <- p_data$sstky - rlag(p_data$sstky,1)
    p_data[p_data$fqtr==1, 'eqissq']  <- p_data[p_data$fqtr==1, 'sstky']
    p_data$prstkcq <- p_data$prstkcy - rlag(p_data$prstkcy,1)
    p_data[p_data$fqtr==1, 'prstkcq'] <- p_data[p_data$fqtr==1, 'prstkcy']

    # Net issuances
    p_data$neqissq <- p_data$eqissq - p_data$prstkcq
    p_data$neqissq_real <- 100*p_data$neqissq / p_data$GVADef
    # Compute ratios to assets (right away in percentages)
    p_data$neqat_rat_realh <- 100*p_data$neqissq_real / rlag(p_data$atq_real,1) # Issuance/lagged assets ratio
}

# Compute quarterly and annual sales growth
p_data$dlsaleq_real  <- p_data$lsaleq_real - rlag(p_data$lsaleq_real,1)
p_data$dylsaleq_real <- p_data$lsaleq_real - rlag(p_data$lsaleq_real,4)
# Compute kinv_rat (based on perpetual inventory method k, right away in percentages)
p_data$kinv_rat_real    <- 100*p_data$capxq_real/rlag(p_data$k_stock, 1)
# Compute "real" inv_rat (right away in percentages)
p_data$inv_rat_realh   <- 100*p_data$capxq_real/rlag(p_data$ppentq_real,1)
# Create investment rate based on difference of CAPXQ and SPPEQ
p_data$sinv_rat_realh   <- 100*(p_data$capxq_real-p_data$sppeq_real)/rlag(p_data$ppentq_real,1)
p_data[is.na(p_data$sinv_rat_realh) & !is.na(p_data$inv_rat_realh),"sinv_rat_realh"]   <- p_data[is.na(p_data$sinv_rat_realh) & !is.na(p_data$inv_rat_realh),"inv_rat_realh"] # If capital sales SPPEQ missing, use gross (CAPXQ) measure

# Date variables into format
p_data$cqtr_num <- as.double(as.character(p_data$cqtr_num))
p_data$datacqtr <- as.factor(p_data$datacqtr)
# Compute age
p_data['age'] = p_data['cqtr_num'] - p_data['incdate_num']

########
# Preparation of panel for regressions
########
library(lfe)
options(lfe.eps=1e-7) # Set tolerance in lfe package
# Copy data to modify
p_data_capx_q <- p_data
# Choose SHOCK MEASURE as JK's FF4 and RESCALE with sd
p_data_capx_q$sffr <- p_data_capx_q$sffr_ff4 / sd(aggregate(sffr_ff4 ~ cqtr_num, data[(data$cqtr_num>=1990) & (data$cqtr_num<=2017),], function(x) mean(x))[,'sffr_ff4'])
# Recount observations and only keep firms with more than 40 post-1990
p_data_capx_q <- p_data_capx_q[p_data_capx_q$cqtr_num >= 1990.0,]
p_data_capx_q$obsc <- rep(count(p_data_capx_q,'gvkey')$freq,count(p_data_capx_q,'gvkey')$freq)
p_data_capx_q <- p_data_capx_q[p_data_capx_q$obsc >= 40,]

# Winsorize ratios
win_val_q = 0.01
if (TRUE){
    # Two-sided
    winsor_keep_fn_q <- function(x) x >= quantile(x,win_val_q,na.rm=TRUE) & x <= quantile(x,1.0-win_val_q,na.rm=TRUE)
    #
    p_data_capx_q <- ddply(p_data_capx_q, c("cqtr_num"), transform, keep_cheat_rat_dm = winsor_keep_fn_q(dlsaleq_real))
    p_data_capx_q[!p_data_capx_q$keep_cheat_rat_dm | is.na(p_data_capx_q$keep_cheat_rat_dm) ,c("dlsaleq_real")] <- NA
    p_data_capx_q <- ddply(p_data_capx_q, c("cqtr_num"), transform, keep_cat_rat = winsor_keep_fn_q(kinv_rat_real))
    p_data_capx_q[!p_data_capx_q$keep_cat_rat | is.na(p_data_capx_q$keep_cat_rat) ,c("kinv_rat_real")] <- NA
    p_data_capx_q <- ddply(p_data_capx_q, c("cqtr_num"), transform, keep_cat_rat = winsor_keep_fn_q(inv_rat_realh))
    p_data_capx_q[!p_data_capx_q$keep_cat_rat | is.na(p_data_capx_q$keep_cat_rat) ,c("inv_rat_realh")] <- NA
    p_data_capx_q <- ddply(p_data_capx_q, c("cqtr_num"), transform, keep_cat_rat = winsor_keep_fn_q(neqat_rat_realh))
    p_data_capx_q[!p_data_capx_q$keep_cat_rat | is.na(p_data_capx_q$keep_cat_rat) ,c("neqat_rat_realh")] <- NA
    p_data_capx_q <- ddply(p_data_capx_q, c("cqtr_num"), transform, keep_cat_rat = winsor_keep_fn_q(sinv_rat_realh))
    p_data_capx_q[!p_data_capx_q$keep_cat_rat | is.na(p_data_capx_q$keep_cat_rat) ,c("sinv_rat_realh")] <- NA

    # One-sided (top)
    winsor_keep_fn_q <- function(x) x <= quantile(x,1.0-win_val_q,na.rm=TRUE)
    p_data_capx_q <- ddply(p_data_capx_q, c("cqtr_num"), transform, keep_lev_rat = winsor_keep_fn_q(lev_rat))
    p_data_capx_q[!p_data_capx_q$keep_lev_rat | is.na(p_data_capx_q$keep_lev_rat) ,c("lev_rat")] <- NA
    p_data_capx_q <- ddply(p_data_capx_q, c("cqtr_num"), transform, keep_cheat_rat = winsor_keep_fn_q(cheat_rat))
    p_data_capx_q[!p_data_capx_q$keep_cheat_rat | is.na(p_data_capx_q$keep_cheat_rat) ,c("cheat_rat")] <- NA
    p_data_capx_q <- ddply(p_data_capx_q, c("cqtr_num"), transform, keep_cat_rat = winsor_keep_fn_q(tobin_q))
    p_data_capx_q[!p_data_capx_q$keep_cat_rat | is.na(p_data_capx_q$keep_cat_rat) ,c("tobin_q")] <- NA
    p_data_capx_q <- ddply(p_data_capx_q, c("cqtr_num"), transform, keep_cat_rat = winsor_keep_fn_q(mtover))
    p_data_capx_q[!p_data_capx_q$keep_cat_rat | is.na(p_data_capx_q$keep_cat_rat) ,c("mtover")] <- NA
    p_data_capx_q <- ddply(p_data_capx_q, c("cqtr_num"), transform, keep_cat_rat = winsor_keep_fn_q(ltat_rat))
    p_data_capx_q[!p_data_capx_q$keep_cat_rat | is.na(p_data_capx_q$keep_cat_rat) ,c("ltat_rat")] <- NA
    p_data_capx_q <- ddply(p_data_capx_q, c("cqtr_num"), transform, keep_cat_rat = winsor_keep_fn_q(olqat_rat))
    p_data_capx_q[!p_data_capx_q$keep_cat_rat | is.na(p_data_capx_q$keep_cat_rat) ,c("olqat_rat")] <- NA
    p_data_capx_q <- ddply(p_data_capx_q, c("cqtr_num"), transform, keep_cat_rat = winsor_keep_fn_q(ppentqat_rat))
    p_data_capx_q[!p_data_capx_q$keep_cat_rat | is.na(p_data_capx_q$keep_cat_rat) ,c("ppentqat_rat")] <- NA

    # Reorder data
    p_data_capx_q <- p_data_capx_q[order(p_data_capx_q$gvkey, p_data_capx_q$cqtr_num), ]
    p_data_capx_q <- pdata.frame(p_data_capx_q, index=c("gvkey", "cqtr_num"))
}

# Compute average standard deviation of turnover across time and normalize turnover
if (TRUE){
    t_meanmtover = aggregate(mtover ~ cqtr_num, p_data_capx_q, function(x) mean(x))
    t_sdmtover = aggregate(mtover ~ cqtr_num, p_data_capx_q, function(x) sd(x))
    avg_sdmtover = mean(t_sdmtover[(as.double(as.character(t_sdmtover$cqtr_num))>= 1990.0) & (as.double(as.character(t_sdmtover$cqtr_num)) <= 2017.00), 'mtover'])
    # Normalize turnover
    p_data_capx_q$mtover <- p_data_capx_q$mtover / avg_sdmtover
}

# Create log(q) and "stobin_q" which is simply q, but dropping all q values above 10 -- e.g. as in Gomes (2001)
p_data_capx_q$ltobin_q <- 100*log(p_data_capx_q$tobin_q)
p_data_capx_q["stobin_q"] <- p_data_capx_q$tobin_q
p_data_capx_q[(p_data_capx_q$stobin_q > 10) | is.na(p_data_capx_q$stobin_q),"stobin_q"] <- NA
# Create log investment rates
p_data_capx_q$linv_rat_realh <- 100*log(p_data_capx_q$inv_rat_realh)
p_data_capx_q[p_data_capx_q$linv_rat_realh==-Inf | is.na(p_data_capx_q$linv_rat_realh),"linv_rat_realh"] <- NA

# For balance sheet (log)-level responses, simply drop the bottom end 0.1%, and any -Inf values
if (TRUE){
    w_cut_balsh = 0.001
    p_data_capx_q[(p_data_capx_q$latq < quantile(p_data_capx_q$latq, w_cut_balsh, na.rm=TRUE)) | is.na(p_data_capx_q$latq),"latq"] <- NA
    p_data_capx_q[(p_data_capx_q$lltq < quantile(p_data_capx_q$lltq, w_cut_balsh, na.rm=TRUE)) | is.na(p_data_capx_q$lltq),"lltq"] <- NA
    p_data_capx_q[(p_data_capx_q$lppentq < quantile(p_data_capx_q$lppentq, w_cut_balsh, na.rm=TRUE)) | is.na(p_data_capx_q$lppentq),"lppentq"] <- NA
    p_data_capx_q[(p_data_capx_q$ldttq < quantile(p_data_capx_q$ldttq, w_cut_balsh, na.rm=TRUE)) | is.na(p_data_capx_q$ldttq),"ldttq"] <- NA
    p_data_capx_q[p_data_capx_q$ldttq==-Inf | is.na(p_data_capx_q$ldttq),"ldttq"] <- NA
}

# Split HIGH/LOW cash
if (TRUE){
    # Demean at time level
    p_data_capx_q <- transform(p_data_capx_q, cheat_rat_tdm = cheat_rat-ave(cheat_rat, datacqtr, FUN=function(x) median(x, na.rm=T)))
    # With indicator
    p_data_capx_q["lliq_ind"] <- NA
    p_data_capx_q[p_data_capx_q$cheat_rat_tdm >= 0 | is.na(p_data_capx_q$cheat_rat_tdm), "lliq_ind"] <- 0
    p_data_capx_q[p_data_capx_q$cheat_rat_tdm < 0 | is.na(p_data_capx_q$cheat_rat_tdm), "lliq_ind"] <- 1
    p_data_capx_q[is.na(p_data_capx_q$cheat_rat_tdm), "lliq_ind"] <- NA

    p_data_capx_q <- p_data_capx_q[order(p_data_capx_q$gvkey, p_data_capx_q$cqtr_num), ]
    p_data_capx_q <- pdata.frame(p_data_capx_q, index=c("gvkey", "cqtr_num"))
    p_data_capx_q$cqtr_num <- as.double(as.character(p_data_capx_q$cqtr_num))
}
# Split HIGH/LOW turnover
if (TRUE){
    # Demean at time level
    p_data_capx_q <- transform(p_data_capx_q, mtover_tdm = mtover-ave(mtover, datacqtr, FUN=function(x) median(x, na.rm=T)))
    # With indicator
    p_data_capx_q["Hmtover_ind"] <- NA
    p_data_capx_q[p_data_capx_q$mtover_tdm >= 0 | is.na(p_data_capx_q$mtover_tdm), "Hmtover_ind"] <- 1
    p_data_capx_q[p_data_capx_q$mtover_tdm < 0 | is.na(p_data_capx_q$mtover_tdm), "Hmtover_ind"] <- 0
    p_data_capx_q[is.na(p_data_capx_q$mtover_tdm), "Hmtover_ind"] <- NA

    p_data_capx_q <- p_data_capx_q[order(p_data_capx_q$gvkey, p_data_capx_q$cqtr_num), ]
    p_data_capx_q <- pdata.frame(p_data_capx_q, index=c("gvkey", "cqtr_num"))
    p_data_capx_q$cqtr_num <- as.double(as.character(p_data_capx_q$cqtr_num))
}
# Create (gvkey X lliq)-dummy and (sic2_datacqtr X lliq)-dummy
if (TRUE){
    # gvkey
    p_data_capx_q$gvkeyXlliq <- paste(as.character(p_data_capx_q$gvkey), as.character(p_data_capx_q$lliq_ind), sep=",")
    p_data_capx_q$gvkeyXlliq <- as.factor(p_data_capx_q$gvkeyXlliq)
    # sic2_datacqtr
    p_data_capx_q$sic2_datacqtrXlliq <- paste(as.character(p_data_capx_q$sic2_datacqtr), as.character(p_data_capx_q$lliq_ind), sep=",")
    p_data_capx_q$sic2_datacqtrXlliq <- as.factor(p_data_capx_q$sic2_datacqtrXlliq)
}

###
# Reduced-form OLS regression (Figure 1)
###
# Turn warnings off (=-1) or on (=0)
options(warn=-1)

# START dep_var loop
for (dep_var in c("ltobin_q", "neqat_rat_realh", "linv_rat_realh")){
    for (selector_q_loop in c("mtover")){
# Initialize
selector_q  = selector_q_loop

# Specification parameters
ylim_man_ind = FALSE
ylim_man     = c(-8.5,2.1) # Manual limits on y-axis
str_end = ""
ALLin_ind   = FALSE # No other cross-term controls
fe_ind      = TRUE   # Firm fixed effect
indt_ind    = TRUE # Industry x time indicator
ind_ind     = FALSE       # Just industry indicator
fc_ind      = TRUE   # Firm controls
cluftn_ind  = TRUE # Cluster standard errors at (timeXind,firm)-level
cluft_ind   = TRUE # Cluster standard errors at (time,firm)-level
clu_ind     = TRUE # Cluster standard errors at time-level
wd_ind      = FALSE  # Shocks weighted depending on day in quarter
ci_lev_set = 0.95

# Lag extra back
lplus_capx_q = 0
# Loop
lag_values <- seq(0,12,1)
# Create collectors
betas_vec_q <- matrix(0,length(lag_values),1)
ci_vec_q <- array(0,dim=c(length(lag_values),2))

# Add correct string at end of filenames
if (fe_ind) {str_end <- paste(str_end,"_FE", sep="")}
if (indt_ind) {str_end <- paste(str_end,"_INDT", sep="")} else{if (ind_ind) {str_end <- paste(str_end,"_IND", sep="")} }
if (fc_ind) {str_end <- paste(str_end,"_FC", sep="")}
if (cluftn_ind) {str_end <- paste(str_end,"_CLUftn", sep="")} else { if (cluft_ind) {str_end <- paste(str_end,"_CLUft", sep="")} else {if (clu_ind) {str_end <- paste(str_end,"_CLU", sep="")}} }
if (wd_ind) {str_end <- paste(str_end,"_wd", sep=""); shock_name_ext<-"_wd"} else {str_end <- paste(str_end,"_nwd", sep=""); shock_name_ext<-""}
if (ALLin_ind) {str_end <- paste(str_end,"_ALLin", sep=""); allmodel_list <- list()}

### Start loop
ctr <- 0
for (lplus_capx_q in lag_values){
ctr <- ctr +1
print('#####')
print(lplus_capx_q)
print('#####')
lag0_capx_q=0 + lplus_capx_q # first lag of sffr
lag1_capx_q=0 + lplus_capx_q # last lag of sffr
lagc0_capx_q=0 + lplus_capx_q # lag of earlier controls
lagc1_capx_q=1 + lplus_capx_q # lag of later controls

# Extra firm controls to include interacted, if requested
firm_contrs = c("rlag(cheat_rat, lagc1_capx_q)", "rlag(lev_rat, lagc1_capx_q)", "rlag(latq, lagc1_capx_q)")

### Cross-terms with monetary shock
if (ALLin_ind){
    controls_ct_capx_q =c(paste("rlag(",selector_q,", lagc1_capx_q + 0)", sep=""), firm_contrs)
} else {
    controls_ct_capx_q =c(paste("rlag(",selector_q,", lagc1_capx_q + 0)", sep=""))
}

### Non-cross-terms
if (ALLin_ind){
    controls_nct_capx_q = c(paste("rlag(",selector_q,", lagc1_capx_q + 0)", sep=""), paste("rlag(",dep_var,", lagc1_capx_q:(lagc1_capx_q))", sep=""), firm_contrs)
} else {
    if (!fc_ind) {
        # No firm controls
        controls_nct_capx_q =c(paste("rlag(",selector_q,", lagc1_capx_q + 0)", sep=""))
    } else {
        # With firm controls
        controls_nct_capx_q =c(paste("rlag(",selector_q,", lagc1_capx_q + 0)", sep=""), paste("rlag(",dep_var,", lagc1_capx_q:(lagc1_capx_q))", sep=""), "rlag(latq, lagc1_capx_q)", "rlag(cheat_rat, lagc1_capx_q)", "rlag(lev_rat, lagc1_capx_q)")
    }
}

# Set up LHS variable and panel for running regression
p_data_capx_q["diff_dep_var"] <- p_data_capx_q[dep_var]
p_data_capx_q_used <- p_data_capx_q

### Control dummies
# If requested, do industry*time fixed effects instead
if (indt_ind) {
    dummies_capx_q = c("sic2_datacqtr")
} else {
    if (ind_ind) {
        # Or if not, and only industry dummies requested
        dummies_capx_q = c("datacqtr", "sic1")
    } else {
        # Otherwise no industry dummies
        dummies_capx_q = c("datacqtr")
    }
}
# If requested, add firm fixed effects
if (fe_ind) {
    dummies_capx_q = c(dummies_capx_q, "gvkey")
}

# Formula
if (cluftn_ind) {
    formla_capx_q <- formula(paste("diff_dep_var ~ rlag(sffr",shock_name_ext,", lag0_capx_q:lag1_capx_q):(", paste(controls_ct_capx_q, collapse=" + "), ") + ", paste(controls_nct_capx_q, collapse=" + "), "  | ", paste(dummies_capx_q, collapse=" + "), "|0 | sic3_datacqtr + gvkey", sep="" ))
} else {
    if (cluft_ind) {
        formla_capx_q <- formula(paste("diff_dep_var ~ rlag(sffr",shock_name_ext,", lag0_capx_q:lag1_capx_q):(", paste(controls_ct_capx_q, collapse=" + "), ") + ", paste(controls_nct_capx_q, collapse=" + "), "  | ", paste(dummies_capx_q, collapse=" + "), "|0 | datacqtr + gvkey", sep="" ))
    } else {
        if (clu_ind){
            formla_capx_q <- formula(paste("diff_dep_var ~ rlag(sffr",shock_name_ext,", lag0_capx_q:lag1_capx_q):(", paste(controls_ct_capx_q, collapse=" + "), ") + ", paste(controls_nct_capx_q, collapse=" + "), "  | ", paste(dummies_capx_q, collapse=" + "), "|0 | datacqtr", sep="" ))
        } else{
            formla_capx_q <- formula(paste("diff_dep_var ~ rlag(sffr",shock_name_ext,", lag0_capx_q:lag1_capx_q):(", paste(controls_ct_capx_q, collapse=" + "), ") + ", paste(controls_nct_capx_q, collapse=" + "), "  | ", paste(dummies_capx_q, collapse=" + "), "|0 | gvkey" , sep=""))
        }
    }
}

# Estimate
reg_capx_q <- felm(formla_capx_q, data=p_data_capx_q_used, na.action=na.exclude)
reg_sum_capx_q <- summary(reg_capx_q)
message(print(reg_sum_capx_q))

reg_coefs_capx_q <- coef(reg_sum_capx_q)
# Save coefs
betas_vec_q[ctr,1] <- reg_coefs_capx_q[, "Estimate"][c(paste("rlag(sffr",shock_name_ext,", lag0_capx_q:lag1_capx_q):rlag(",selector_q,", lagc1_capx_q + 0)", sep=""))]
ci_vec_q[ctr,] <- confint(reg_capx_q, level=ci_lev_set)[c(paste("rlag(sffr",shock_name_ext,", lag0_capx_q:lag1_capx_q):rlag(",selector_q,", lagc1_capx_q + 0)", sep="")),]

} # End of estimations loop

# Create folder for figures
if (!dir.exists("./saved_figs")) {dir.create("./saved_figs")}
if (!dir.exists("./saved_figs/_fig1")) {dir.create("./saved_figs/_fig1")}
# Plot results for main variable
pdf(paste("./saved_figs/_fig1/slevs_betas_vec_q(",dep_var,")(",selector_q,")",str_end,".pdf",sep=""), width=8, height=7)
par(mar = c(5,5.3,2,2))
cl<- rainbow(1, start=0.7)
cl <- c("black")
markers <- rep(c(4, 19), 1)
if (ylim_man_ind){
    plot(lag_values, betas_vec_q[,1], las=1, type="n", ylim=ylim_man, xlab = "Quarters (h)", ylab=expression(gamma["h"]^x), cex.lab=font_lab_val, cex.axis=font_ax_val)
} else {
    plot(lag_values, betas_vec_q[,1], type="n", ylim=c(-0.0002+min(0.00, min(ci_vec_q[,1])) , 0.0002+max(0,max(ci_vec_q[,2]))), xlab = "Quarters (h)", ylab=expression(gamma["h"]), cex.lab=font_lab_val, cex.axis=font_ax_val)
}
lines(lag_values,betas_vec_q[,1], type="o", col=cl[1], pch=markers[1], lwd=linew_main); lines(lag_values, ci_vec_q[,1], lty=2, col=cl[1], lwd=linew_ci); lines(lag_values,ci_vec_q[,2], lty=2, col=cl[1], lwd=linew_ci);
abline(h=0.0); grid (NULL,NULL, lty = 6, col = "cornsilk2")
dev.off()

} # End of "selector_q" loop
} # End of "dep_var" loop


########
# Reduced-form OLS regression for high/low liquidity ratios (Figure 2)
########

# Iterate over figure numbers
fig_list = c("_fig2", "_fig13")

for (fig_name in fig_list){
    # Set dependent variables to be considered in Figure
    if (fig_name %in% c("_fig2")){
        dep_var_list = c("ltobin_q","neqat_rat_realh","linv_rat_realh")
    } else if (fig_name %in% c("_fig13")) {
        dep_var_list = c("ltobin_q")
    }

    # Figure 13 uses JKpm shocks
    if (fig_name %in% c("_fig13")){
        p_data_capx_q$sffr <- p_data_capx_q$sffr_JKpmt / sd(aggregate(sffr_JKpmt ~ cqtr_num, data[(data$cqtr_num>=1990) & (data$cqtr_num<=2017),], function(x) mean(x))[,'sffr_JKpmt'])
    }

# START dep_var loop
for (dep_var in dep_var_list){
    for (selector_q_loop in c("mtover")){
# Initialize
selector_q  = selector_q_loop

# Specification parameters
ylim_man_ind = FALSE
ylim_man     = c(-10,2.1)  # Manual limits on y-axis
str_end = ""
ALLin_ind   = FALSE # No other cross-term controls
fe_ind      = TRUE   # Firm fixed effect
indt_ind    = TRUE # Industry x time indicator
ind_ind     = FALSE       # Just industry indicator
fc_ind      = TRUE   # Firm controls
cluftn_ind  = TRUE # Cluster standard errors at (timeXind,firm)-level
cluft_ind   = TRUE # Cluster standard errors at (time,firm)-level
clu_ind     = TRUE # Cluster standard errors at time-level
wd_ind      = FALSE  # Weighted shocks or not
ci_true_ind = TRUE # Compute correct confidence intervals for sums of parameters
wd_ind      = FALSE  # Shocks weighted depending on day in quarter
ci_lev_set = 0.95
# Labels on figures
hliq_label   = "High liq"
lliq_label = "Low liq"

# Lag extra back
lplus_capx_q = 0
# Loop
lag_values <- seq(0,12,1)
# Create collectors
betas_vec_q <- matrix(0,length(lag_values),1)
ci_vec_q <- array(0,dim=c(length(lag_values),2))
yb_vec_q <- matrix(0,length(lag_values),1)
yb_ci_vec_q <- array(0,dim=c(length(lag_values),2))
ylevb_vec_q <- matrix(0,length(lag_values),1)
ylevb_ci_vec_q <- array(0,dim=c(length(lag_values),2))

# Add correct string at end of filenames
if (fe_ind) {str_end <- paste(str_end,"_FE", sep="")}
if (indt_ind) {str_end <- paste(str_end,"_INDT", sep="")} else{if (ind_ind) {str_end <- paste(str_end,"_IND", sep="")} }
if (fc_ind) {str_end <- paste(str_end,"_FC", sep="")}
if (cluftn_ind) {str_end <- paste(str_end,"_CLUftn", sep="")} else { if (cluft_ind) {str_end <- paste(str_end,"_CLUft", sep="")} else {if (clu_ind) {str_end <- paste(str_end,"_CLU", sep="")}} }
if (wd_ind) {str_end <- paste(str_end,"_wd", sep=""); shock_name_ext<-"_wd"} else {str_end <- paste(str_end,"_nwd", sep=""); shock_name_ext<-""}
if (ALLin_ind) {str_end <- paste(str_end,"_ALLin", sep=""); allmodel_list <- list()}

### Start loop
ctr <- 0
for (lplus_capx_q in lag_values){
ctr <- ctr +1
print('#####')
print(lplus_capx_q)
print('#####')
lag0_capx_q=0 + lplus_capx_q # first lag of dffr
lag1_capx_q=0 + lplus_capx_q # last lag of dffr
lagc0_capx_q=0 + lplus_capx_q # lag of earlier controls
lagc1_capx_q=1 + lplus_capx_q # lag of later controls

# Extra firm controls to include interacted, if requested
firm_contrs = c("rlag(cheat_rat, lagc1_capx_q)", "rlag(lev_rat, lagc1_capx_q)", "rlag(latq, lagc1_capx_q)")

### Cross-terms with monetary shock
if (ALLin_ind){
    controls_ct_capx_q =c(paste("rlag(",selector_q,", lagc1_capx_q + 0)", sep=""), "rlag(lliq_ind, lagc1_capx_q)", paste("rlag(",selector_q,", lagc1_capx_q + 0):rlag(lliq_ind, lagc1_capx_q)", sep=""), firm_contrs, paste("rlag(lliq_ind, lagc1_capx_q):(", paste(firm_contrs, collapse=" + "), ")", sep=""))
} else {
    controls_ct_capx_q =c(paste("rlag(",selector_q,", lagc1_capx_q + 0)", sep=""), "rlag(lliq_ind, lagc1_capx_q)", paste("rlag(",selector_q,", lagc1_capx_q + 0):rlag(lliq_ind, lagc1_capx_q)", sep=""))
}

### Non-cross-terms
if (ALLin_ind){
    controls_nct_capx_q = c(paste("rlag(",selector_q,", lagc1_capx_q + 0)", sep=""), "rlag(lliq_ind, lagc1_capx_q)", paste("rlag(",selector_q,", lagc1_capx_q + 0):rlag(lliq_ind, lagc1_capx_q)", sep=""), firm_contrs, paste("rlag(lliq_ind, lagc1_capx_q):(", paste(firm_contrs, collapse=" + "), ")", sep=""), paste("rlag(",dep_var,", lagc1_capx_q:(lagc1_capx_q+0))", sep=""), paste("rlag(lliq_ind, lagc1_capx_q):rlag(",dep_var,", lagc1_capx_q:(lagc1_capx_q+0))", sep=""))
} else {
    # With firm controls
    controls_nct_capx_q =c(paste("rlag(",selector_q,", lagc1_capx_q + 0)", sep=""), "rlag(lliq_ind, lagc1_capx_q)", paste("rlag(",selector_q,", lagc1_capx_q + 0):rlag(lliq_ind, lagc1_capx_q)", sep=""), paste("rlag(",dep_var,", lagc1_capx_q:(lagc1_capx_q+0))", sep=""), paste("rlag(lliq_ind, lagc1_capx_q):rlag(",dep_var,", lagc1_capx_q:(lagc1_capx_q+0))", sep=""), "rlag(lev_rat, lagc1_capx_q)", "rlag(latq, lagc1_capx_q)", "rlag(cheat_rat, lagc1_capx_q)", "rlag(lliq_ind, lagc1_capx_q):rlag(lev_rat, lagc1_capx_q)", "rlag(lliq_ind, lagc1_capx_q):rlag(latq, lagc1_capx_q)", "rlag(lliq_ind, lagc1_capx_q):rlag(cheat_rat, lagc1_capx_q)")
}

# Set up LHS variable and panel for running regression

# Create differences
p_data_capx_q["diff_dep_var"] <- p_data_capx_q[dep_var]
p_data_capx_q_used <- p_data_capx_q

### Control dummies
# Create correctly lagged (sic2_datacqtr X lliq) FE
p_data_capx_q_used$sic2_datacqtrXlliq_cur <- rlag(p_data_capx_q_used$sic2_datacqtrXlliq,lagc1_capx_q) # SIC2 X time X lliq FE
if (indt_ind) {
    dummies_capx_q = c("sic2_datacqtrXlliq_cur")
} else {
    if (ind_ind) {
        # Or if not, and only industry dummies requested
        dummies_capx_q = c("datacqtr", "sic1")
    } else {
        # Otherwise no industry dummies
        dummies_capx_q = c("datacqtr")
    }
}
# Create correctly lagged (firm X lliq) FE
p_data_capx_q_used$gvkeyXlliq_cur <- rlag(p_data_capx_q_used$gvkeyXlliq,lagc1_capx_q) # Firm X lliq FE
if (fe_ind) {
    dummies_capx_q = c(dummies_capx_q, "gvkeyXlliq_cur")
}

# Formula
if (cluftn_ind) {
    formla_capx_q <- formula(paste("diff_dep_var ~ rlag(sffr",shock_name_ext,", lag0_capx_q:lag1_capx_q):(", paste(controls_ct_capx_q, collapse=" + "), ") + ", paste(controls_nct_capx_q, collapse=" + "), "  | ", paste(dummies_capx_q, collapse=" + "), "|0 | sic3_datacqtr + gvkey", sep="" ))
} else {
    if (cluft_ind) {
        formla_capx_q <- formula(paste("diff_dep_var ~ rlag(sffr",shock_name_ext,", lag0_capx_q:lag1_capx_q):(", paste(controls_ct_capx_q, collapse=" + "), ") + ", paste(controls_nct_capx_q, collapse=" + "), "  | ", paste(dummies_capx_q, collapse=" + "), "|0 | datacqtr + gvkey", sep="" ))
    } else {
        if (clu_ind){
            formla_capx_q <- formula(paste("diff_dep_var ~ rlag(sffr",shock_name_ext,", lag0_capx_q:lag1_capx_q):(", paste(controls_ct_capx_q, collapse=" + "), ") + ", paste(controls_nct_capx_q, collapse=" + "), "  | ", paste(dummies_capx_q, collapse=" + "), "|0 | datacqtr", sep="" ))
        } else{
            formla_capx_q <- formula(paste("diff_dep_var ~ rlag(sffr",shock_name_ext,", lag0_capx_q:lag1_capx_q):(", paste(controls_ct_capx_q, collapse=" + "), ") + ", paste(controls_nct_capx_q, collapse=" + "), "  | ", paste(dummies_capx_q, collapse=" + "), "|0 | gvkey" , sep=""))
        }
    }
}

# Estimate
reg_capx_q <- felm(formla_capx_q, data=p_data_capx_q_used, na.action=na.exclude)
reg_sum_capx_q <- summary(reg_capx_q)
message(print(reg_sum_capx_q))

reg_coefs_capx_q <- coef(reg_sum_capx_q)
# Save coefs
betas_vec_q[ctr,1] <- reg_coefs_capx_q[, "Estimate"][c(paste("rlag(sffr",shock_name_ext,", lag0_capx_q:lag1_capx_q):rlag(",selector_q,", lagc1_capx_q + 0)", sep=""))]
ci_vec_q[ctr,] <- confint(reg_capx_q, level=ci_lev_set)[c(paste("rlag(sffr",shock_name_ext,", lag0_capx_q:lag1_capx_q):rlag(",selector_q,", lagc1_capx_q + 0)", sep="")),]
# Save coefs for lliq_ind
yb_vec_q[ctr,1] <- reg_coefs_capx_q[, "Estimate"][c(paste("rlag(sffr",shock_name_ext,", lag0_capx_q:lag1_capx_q):rlag(",selector_q,", lagc1_capx_q + 0):rlag(lliq_ind, lagc1_capx_q)", sep=""))]
yb_ci_vec_q[ctr,] <- confint(reg_capx_q, level=ci_lev_set)[c(paste("rlag(sffr",shock_name_ext,", lag0_capx_q:lag1_capx_q):rlag(",selector_q,", lagc1_capx_q + 0):rlag(lliq_ind, lagc1_capx_q)", sep="")),]

# Compute and collect correct ("true") confidence intervals of sums of parameters
if (ci_true_ind){
    crit_val = 1.96
    if (ci_lev_set==0.90){crit_val = 1.645}
    VCV_cur = reg_capx_q$clustervcv
    nK = dim(VCV_cur)[1]
    sel_vec <- numeric(nK)
    lev_ind <- match(paste("rlag(sffr",shock_name_ext,", lag0_capx_q:lag1_capx_q):rlag(",selector_q,", lagc1_capx_q + 0)", sep=""), row.names(VCV_cur))
    yb_ind  <- match(paste("rlag(sffr",shock_name_ext,", lag0_capx_q:lag1_capx_q):rlag(",selector_q,", lagc1_capx_q + 0):rlag(lliq_ind, lagc1_capx_q)", sep=""), row.names(VCV_cur))
    sel_vec[c(lev_ind,yb_ind)] <- c(1,1)
    SE_sum <- sqrt(sel_vec %*% VCV_cur %*% sel_vec)
    # Save
    ylevb_vec_q[ctr,1] <- betas_vec_q[ctr] + yb_vec_q[ctr]
    ylevb_ci_vec_q[ctr,] <- ylevb_vec_q[ctr,1] + c(-1,1) * crit_val * SE_sum
}

} # End of estimations loop

# Create folder for figures
if (!dir.exists(paste("./saved_figs/",fig_name,sep=""))) {dir.create(paste("./saved_figs/",fig_name,sep=""))}

# Plot results for main coefs and coefs + lliq_ind -- with correct confidence intervals
if (ci_true_ind){
    pdf(paste("./saved_figs/",fig_name,"/slevs_betas-yb_vec_q_citrue(",dep_var,")(",selector_q,")",str_end,".pdf",sep=""), width=8, height=7)
    par(mar = c(5,5.3,2,2))
    cl<- rainbow(1, start=0.7)
    cl <- c("blue")
    cl2 <- c("red")
    markers <- rep(c(4, 19), 1)
    if (ylim_man_ind){
        plot(lag_values, betas_vec_q[,1], las=1, type="n", ylim=ylim_man, xlab = "Quarters (h)", ylab=expression(paste(gamma["h"],", ", gamma["h"]+tilde(gamma)["h"])), cex.lab=font_lab_val, cex.axis=font_ax_val)
    } else {
        plot(lag_values, betas_vec_q[,1], type="n", ylim=c(-0.0002+min(0.00, min(ci_vec_q[,1], ylevb_ci_vec_q[,1])) , 0.0002+max(0,max(ci_vec_q[,2], ylevb_ci_vec_q[,2]))), xlab = "Quarters (h)", ylab=expression(paste(gamma["h"],", ", gamma["h"]+tilde(gamma)["h"])), cex.lab=font_lab_val, cex.axis=font_ax_val)
    }
    lines(lag_values,betas_vec_q[,1], type="o", col=cl[1], pch=markers[1], lwd=linew_main); lines(lag_values, ci_vec_q[,1], lty=2, col=cl[1], lwd=linew_ci); lines(lag_values,ci_vec_q[,2], lty=2, col=cl[1], lwd=linew_ci);
    lines(lag_values,ylevb_vec_q[,1], type="o", col=cl2[1], pch=markers[2], lwd=linew_main); lines(lag_values, ylevb_ci_vec_q[,1], lty=2, col=cl2[1], lwd=linew_ci); lines(lag_values, ylevb_ci_vec_q[,2], lty=2, col=cl2[1], lwd=linew_ci);
    abline(h=0.0); grid (NULL,NULL, lty = 6, col = "cornsilk2"); legend("bottomleft", lwd=rep(linew_main,2), lty=rep(2,1), col=c("blue","red"), pch=markers, cex=font_leg_val, legend=c(hliq_label,lliq_label))
    dev.off()
}

} # End of "selector_q" loop
} # End of "dep_var" loop

# If ran Figure 13 and regressions are done, set sffr back to baseline
if (fig_name %in% c("_fig13")){
    p_data_capx_q$sffr <- p_data_capx_q$sffr_ff4 / sd(aggregate(sffr_ff4 ~ cqtr_num, data[(data$cqtr_num>=1990) & (data$cqtr_num<=2017),], function(x) mean(x))[,'sffr_ff4'])
}

} # fig_name loop end


########
# Baseline IV specification (Figure 3)
########

for (dep_var in c("neqat_rat_realh","linv_rat_realh")){
    for (selector_q_loop in c("ltobin_q")){

# Firm-level instrument to be interacted with sffr
fi_inst     = "mtover"

# Initialize
selector_q  = selector_q_loop
ALLin_ind   = FALSE # No other cross-term controls
# Specification parameters
ylim_man_ind = FALSE
ylim_man = c(-5.0,10.0) # Manual limits on y-axis
#n_q = 2
str_end = "_IV"
fe_ind      = TRUE  # Firm fixed effect
indt_ind    = TRUE  # Industry x time indicator
ind_ind     = FALSE # Just industry indicator
fc_ind      = TRUE  # Firm controls
cluftn_ind  = TRUE  # Cluster standard errors at (Industry x time, firm)-level
cluft_ind   = TRUE  # Cluster standard errors at (time,firm)-level
clu_ind     = TRUE  # Cluster standard errors at time-level
wd_ind      = FALSE  # Shocks weighted depending on day in quarter
ci_lev_set  = 0.95
h_q_bar     = 0     # Maximum "h_q" at which to consider q

# Lag extra back
lplus_capx_q = 0
# Loop
lag_values <- seq(0,12,1)
# Create collectors
betas_vec_q <- matrix(0,length(lag_values),1)
ci_vec_q <- array(0,dim=c(length(lag_values),2))

# Add correct string at end of filenames
if (fe_ind) {str_end <- paste(str_end,"_FE", sep="")}
if (indt_ind) {str_end <- paste(str_end,"_INDT", sep="")} else{if (ind_ind) {str_end <- paste(str_end,"_IND", sep="")} }
if (fc_ind) {str_end <- paste(str_end,"_FC", sep="")}
if (cluftn_ind) {str_end <- paste(str_end,"_CLUftn", sep="")} else { if (cluft_ind) {str_end <- paste(str_end,"_CLUft", sep="")} else{if (clu_ind) {str_end <- paste(str_end,"_CLU", sep="")}} }
if (wd_ind) {str_end <- paste(str_end,"_wd", sep=""); shock_name_ext<-"_wd"} else {str_end <- paste(str_end,"_nwd", sep=""); shock_name_ext<-""}
if (ALLin_ind) {str_end <- paste(str_end,"_ALLin", sep=""); allmodel_list <- list()}

### Start loop
ctr <- 0
for (lplus_capx_q in lag_values){
ctr <- ctr +1
print('#####')
print(lplus_capx_q)
print('#####')
lag0_capx_q=0 + lplus_capx_q # first lag of sffr
lag1_capx_q=0 + lplus_capx_q # last lag of sffr
lagc0_capx_q=0 + lplus_capx_q # lag of earlier controls
lagc1_capx_q=1 + lplus_capx_q # lag of later controls
# Also, label the lag as "h", as in math notation, and the corresponding h_q
h   <- lplus_capx_q
h_q <- min(h, h_q_bar)

### Cross-terms, for controlling other channels, as requested
controls_ct_capx_q <- c("1")
if (ALLin_ind){
    controls_ct_capx_q <- c("rlag(lev_rat, lagc1_capx_q)", "rlag(cheat_rat, lagc1_capx_q)", "rlag(latq, lagc1_capx_q)")
} else {
    controls_ct_capx_q <- c("1")
}

# Create explicit cross-terms in tover and sffr
inst_controls_ct_capx_q_all = c(fi_inst)

for (varname in inst_controls_ct_capx_q_all){
    if (wd_ind){
        eval(parse(text=paste("p_data_capx_q$",varname,"Xsffr<-rlag(p_data_capx_q$",varname,",",as.character(lagc1_capx_q),")*rlag(p_data_capx_q$sffr_wd,",as.character(lag0_capx_q),")",sep="")))
    } else {
        eval(parse(text=paste("p_data_capx_q$",varname,"Xsffr<-rlag(p_data_capx_q$",varname,",",as.character(lagc1_capx_q),")*rlag(p_data_capx_q$sffr,",as.character(lag0_capx_q),")",sep="")))
    }
}
# Also, create the "selector_q" i.e. tobin_q variable to be instrumented
eval(parse(text=paste("p_data_capx_q$endg_var <-rlag(p_data_capx_q$",selector_q,",",as.character(h-h_q),")",sep="")))

### Non-cross-terms
# Add tover (fi_inst), lagged dependent variable, and other controls as necessary
if (ALLin_ind){
    controls_nct_capx_q = c(paste("rlag(",fi_inst,", lagc1_capx_q)", sep=""), paste("rlag(",dep_var,", lagc1_capx_q)", sep=""), "rlag(age, lagc1_capx_q)", "rlag(lev_rat, lagc1_capx_q)", "rlag(cheat_rat, lagc1_capx_q)", "rlag(latq, lagc1_capx_q)")
} else {
    if (!fc_ind) {
        # No firm controls
        controls_nct_capx_q =c(paste("rlag(",fi_inst,", lagc1_capx_q)", sep=""), paste("rlag(",dep_var,", lagc1_capx_q)", sep=""))
    } else {
        controls_nct_capx_q =c(paste("rlag(",fi_inst,", lagc1_capx_q)", sep=""), paste("rlag(",dep_var,", lagc1_capx_q)", sep=""))
        # With firm controls
        controls_nct_capx_q =c(paste("rlag(",fi_inst,", lagc1_capx_q)", sep=""), paste("rlag(",dep_var,", lagc1_capx_q)", sep=""), "rlag(latq, lagc1_capx_q)", "rlag(cheat_rat, lagc1_capx_q)", "rlag(lev_rat, lagc1_capx_q)")
    }
}

# Set up LHS variable and panel for running regression
p_data_capx_q["diff_dep_var"] <- p_data_capx_q[dep_var]
p_data_capx_q_used <- p_data_capx_q

### Control dummies
# If requested, do industry*time fixed effects instead
if (indt_ind) {
    dummies_capx_q = c("sic2_datacqtr")
} else {
    if (ind_ind) {
        # Or if not, and only industry dummies requested
        dummies_capx_q = c("datacqtr", "sic1")
    } else {
        # Otherwise no industry dummies
        dummies_capx_q = c("datacqtr")
    }
}
# If requested, add firm fixed effects
if (fe_ind) {
    dummies_capx_q = c(dummies_capx_q, "gvkey")
}

# Formula
if (cluftn_ind){
    formla_capx_q <- formula(paste("diff_dep_var ~ ", paste(controls_nct_capx_q, collapse=" + "), " + rlag(sffr",shock_name_ext,", lag0_capx_q)*(", paste(controls_ct_capx_q, collapse=" + "), ")| ", paste(dummies_capx_q, collapse=" + "), "|(endg_var ~", paste(paste(fi_inst, "Xsffr",sep=""), collapse="+"),") | sic3_datacqtr + gvkey" ))
} else{
    if (cluft_ind) {
        formla_capx_q <- formula(paste("diff_dep_var ~ ", paste(controls_nct_capx_q, collapse=" + "), " + rlag(sffr",shock_name_ext,", lag0_capx_q)*(", paste(controls_ct_capx_q, collapse=" + "), ")| ", paste(dummies_capx_q, collapse=" + "), "|(endg_var ~", paste(paste(fi_inst, "Xsffr",sep=""), collapse="+"),") | datacqtr + gvkey" ))
    } else {
        if (clu_ind){
            formla_capx_q <- formula(paste("diff_dep_var ~ ", paste(controls_nct_capx_q, collapse=" + "), "  | ", paste(dummies_capx_q, collapse=" + "), "|(endg_var ~", paste(paste(fi_inst, "Xsffr",sep=""), collapse="+"),") | datacqtr" ))
        } else{
            formla_capx_q <- formula(paste("diff_dep_var ~ ", paste(controls_nct_capx_q, collapse=" + "), "  | ", paste(dummies_capx_q, collapse=" + "), "|(endg_var ~", paste(paste(fi_inst, "Xsffr",sep=""), collapse="+"),")" ))
        }
    }
}

# Estimate
reg_capx_q <- felm(formla_capx_q, data=p_data_capx_q_used, na.action=na.exclude)
reg_sum_capx_q <- summary(reg_capx_q)
message(print(reg_sum_capx_q))

reg_coefs_capx_q <- coef(reg_sum_capx_q)
# Save coefs
betas_vec_q[ctr,1] <- reg_coefs_capx_q[, "Estimate"][paste("`endg_var(fit)`", sep="")]
ci_vec_q[ctr,] <- confint(reg_capx_q, level=ci_lev_set)[c(paste("`endg_var(fit)`", sep="")),]

} # End of estimations loop

# Create folder for figures
if (!dir.exists("./saved_figs/_fig3")) {dir.create("./saved_figs/_fig3")}
# Plot results for main variable
pdf(paste("./saved_figs/_fig3/slevs_betas_vec_q(",dep_var,")(",selector_q,")(i-",fi_inst,")",str_end,".pdf",sep=""), width=8, height=7)
par(mar = c(5,5,2,2))
cl<- rainbow(1, start=0.7)
cl <- c("black")
markers <- rep(c(4, 19), 1)
if (ylim_man_ind){
    plot(lag_values, betas_vec_q[,1], type="n", ylim=ylim_man, xlab = "Quarters (h)", ylab=expression(gamma["h"]), cex.lab=font_lab_val, cex.axis=font_ax_val)
} else {
    plot(lag_values, betas_vec_q[,1], type="n", ylim=c(-0.0002+min(0.00, min(ci_vec_q[,1])) , 0.0002+max(0,max(ci_vec_q[,2]))), xlab = "Quarters (h)", ylab=expression(gamma["h"]), cex.lab=font_lab_val, cex.axis=font_ax_val)
}
lines(lag_values,betas_vec_q[,1], type="o", col=cl[1], pch=markers[1], lwd=linew_main); lines(lag_values, ci_vec_q[,1], lty=2, col=cl[1], lwd=linew_ci); lines(lag_values,ci_vec_q[,2], lty=2, col=cl[1], lwd=linew_ci);
abline(h=0.0); grid (NULL,NULL, lty = 6, col = "cornsilk2")
dev.off()

} # End of "selector_q" loop
} # End of "dep_var" loop


####
# Main IV specification for high/low liquidity ratios (Figures 4, 5, 8, 9, 10, 11, 12, 13, 14)
####

# For robustness (Figure 10), estimate firm-specific regressions of sales growth onto GDP growth to measure "gdp_beta"
# To save time, load in from backup file below
if (FALSE){
    firm_list   <- unique(as.double(as.character(p_data_capx_q$gvkey)))
    Ngdp_cutoff <- 5 # Only estimate for firms with at least 5 observations of quarterly sales growth
    # Create empty variable to contain "gdp_beta"
    icount = 0
    p_data_capx_q$gdp_beta <- NA
    for (iig in firm_list){
        icount = icount+1
        if (icount %% 100 == 0){ message(print(icount)) }
        pdat_temp       = p_data_capx_q[p_data_capx_q$gvkey==iig,]
        N_dlsaleq_temp  <- sum(!is.na(pdat_temp$dlsaleq_real))
        if (N_dlsaleq_temp >= Ngdp_cutoff){
            gdp_reg_temp    <- lm(dlsaleq_real ~ dlgdp_real, data = pdat_temp)
            p_data_capx_q[p_data_capx_q$gvkey==iig,"gdp_beta"] = coef(gdp_reg_temp)[2]
        }
    }
    # Save as backup
    bu_gdp_beta <- p_data_capx_q$gdp_beta # Backup
    write.csv(bu_gdp_beta, paste("./bu_gdp_beta_temp.csv",sep=""))
} else {
    # Load in backup
    data_gdp_beta = read.csv("./bu_gdp_beta_temp.csv")
    p_data_capx_q["gdp_beta"] <- data_gdp_beta[,2]
}

# Drop top and bottom 1% gdp_beta outliers
gdp_beta_h = quantile(p_data_capx_q$gdp_beta,na.rm=TRUE,0.99)
gdp_beta_l = quantile(p_data_capx_q$gdp_beta,na.rm=TRUE,0.01)
p_data_capx_q[is.na(p_data_capx_q$gdp_beta) | (p_data_capx_q$gdp_beta > gdp_beta_h) | (p_data_capx_q$gdp_beta < gdp_beta_l), 'gdp_beta'] <- NA

# Iterate over figure numbers
fig_list = c("_fig4", "_fig5", "_fig8", "_fig9", "_fig10", "_fig11", "_fig12", "_fig13", "_fig14")

for (fig_name in fig_list){
    # Set dependent variables to be considered in Figure
    if (fig_name %in% c("_fig4", "_fig8", "_fig9", "_fig10", "_fig11", "_fig12", "_fig13")){
        dep_var_list = c("neqat_rat_realh", "linv_rat_realh")
    } else if (fig_name %in% c("_fig5")) {
        dep_var_list = c("latq", "lltq", "ltat_rat", "lppentq", "ppentqat_rat", "cheat_rat", "ldttq", "lev_rat", "olqat_rat")
    } else if (fig_name %in% c("_fig14")) {
        dep_var_list = c("neqat_rat_realh", "sinv_rat_realh")
    }

    # Set the endogenous variable to be considered in figure
    if (fig_name %in% c("_fig4", "_fig5", "_fig8", "_fig9", "_fig10", "_fig11", "_fig12", "_fig13")){
        selector_q = "ltobin_q"
    } else if (fig_name %in% c("_fig14")) {
        selector_q = "stobin_q"
    }

    # Indicator for other cross-terms or not
    if (fig_name %in% c("_fig4", "_fig5", "_fig13", "_fig14")){
        ALLin_ind = FALSE
    } else if (fig_name %in% c("_fig8", "_fig9", "_fig10", "_fig11", "_fig12")) {
        ALLin_ind = TRUE
    }

    # List of cross-terms to include
    if (fig_name %in% c("_fig8")){
        firm_contrs = c("rlag(cheat_rat, lagc1_capx_q)", "rlag(lev_rat, lagc1_capx_q)", "rlag(latq, lagc1_capx_q)")
    } else if (fig_name %in% c("_fig9")){
        firm_contrs = c("rlag(cheat_rat, lagc1_capx_q)", "rlag(lev_rat, lagc1_capx_q)", "rlag(latq, lagc1_capx_q)", "rlag(age, lagc1_capx_q)")
    } else if (fig_name %in% c("_fig10")){
        firm_contrs = c("rlag(gdp_beta, lagc1_capx_q)")
    } else if (fig_name %in% c("_fig11")){
        firm_contrs = c("rlag(MktRF_beta, lagc1_capx_q)", "rlag(SMB_beta, lagc1_capx_q)", "rlag(HML_beta, lagc1_capx_q)")
    } else if (fig_name %in% c("_fig12")){
        firm_contrs = c("rlag(lev_rat, lagc1_capx_q)", "rlag(RETXsd, lagc1_capx_q)")
    }

    # Figure 13 uses JKpm shocks
    if (fig_name %in% c("_fig13")){
        p_data_capx_q$sffr <- p_data_capx_q$sffr_JKpmt / sd(aggregate(sffr_JKpmt ~ cqtr_num, data[(data$cqtr_num>=1990) & (data$cqtr_num<=2017),], function(x) mean(x))[,'sffr_JKpmt'])
    }

# START dep_var loop
for (dep_var in dep_var_list){

# Firm-level instrument to be interacted with sffr
fi_inst     = "mtover"

# Specification parameters
ylim_man_ind = FALSE
ylim_man = c(-0.5,1.0) # Manual limits on y-axis
str_end = "_IV"
fe_ind      = TRUE  # Firm fixed effect
indt_ind    = TRUE  # Industry x time indicator
ind_ind     = FALSE # Just industry indicator
fc_ind      = TRUE  # Firm controls
cluftn_ind  = TRUE  # Cluster standard errors at (Industry x time, firm)-level
cluft_ind   = TRUE  # Cluster standard errors at (time,firm)-level
clu_ind     = TRUE  # Cluster standard errors at time-level
wd_ind      = FALSE # Weighted shocks or not
ci_true_ind = TRUE # Compute correct confidence intervals for sums of parameters
wd_ind      = FALSE  # Shocks weighted depending on day in quarter
ci_lev_set  = 0.95
h_q_bar     = 0     # Maximum "h_q" at which to consider q
# Labels on figures
hliq_label   = "High liq"
lliq_label = "Low liq"

# Lag extra back
lplus_capx_q = 0
# Loop
lag_values <- seq(0,12,1)
# Create collectors
betas_vec_q <- matrix(0,length(lag_values),1)
ci_vec_q <- array(0,dim=c(length(lag_values),2))
yb_vec_q <- matrix(0,length(lag_values),1)
yb_ci_vec_q <- array(0,dim=c(length(lag_values),2))
ylevb_vec_q <- matrix(0,length(lag_values),1)
ylevb_ci_vec_q <- array(0,dim=c(length(lag_values),2))

# Add correct string at end of filenames
if (fe_ind) {str_end <- paste(str_end,"_FE", sep="")}
if (indt_ind) {str_end <- paste(str_end,"_INDT", sep="")} else{if (ind_ind) {str_end <- paste(str_end,"_IND", sep="")} }
if (fc_ind) {str_end <- paste(str_end,"_FC", sep="")}
if (cluftn_ind) {str_end <- paste(str_end,"_CLUftn", sep="")} else { if (cluft_ind) {str_end <- paste(str_end,"_CLUft", sep="")} else{if (clu_ind) {str_end <- paste(str_end,"_CLU", sep="")}} }
if (wd_ind) {str_end <- paste(str_end,"_wd", sep=""); shock_name_ext<-"_wd"} else {str_end <- paste(str_end,"_nwd", sep=""); shock_name_ext<-""}
if (ALLin_ind) {str_end <- paste(str_end,"_ALLin", sep=""); allmodel_list <- list()}

### Start loop
ctr <- 0
for (lplus_capx_q in lag_values){
ctr <- ctr +1
print('#####')
print(lplus_capx_q)
print('#####')
lag0_capx_q=0 + lplus_capx_q # first lag of sffr
lag1_capx_q=0 + lplus_capx_q # last lag of sffr
lagc0_capx_q=0 + lplus_capx_q # lag of earlier controls (cash flows?)
lagc1_capx_q=1 + lplus_capx_q # lag of later controls
# Also, label the lag as "h", like in math notation, and the corresponding h_q
h   <- lplus_capx_q
h_q <- min(h, h_q_bar)

if (ALLin_ind){
    controls_ct_capx_q <- c("rlag(lliq_ind, lagc1_capx_q)", firm_contrs, paste("rlag(lliq_ind, lagc1_capx_q):(", paste(firm_contrs, collapse=" + "), ")", sep=""))
} else {
    controls_ct_capx_q <- c("rlag(lliq_ind, lagc1_capx_q)")
}

# Create explicit cross-terms in (tover,sffr), and (tover,sffr,lliq_ind)
inst_controls_ct_capx_q_all = c(fi_inst)

for (varname in inst_controls_ct_capx_q_all){
    if (wd_ind){
        eval(parse(text=paste("p_data_capx_q$",varname,"Xsffr<-rlag(p_data_capx_q$",varname,",",as.character(lagc1_capx_q),")*rlag(p_data_capx_q$sffr_wd,",as.character(lag0_capx_q),")",sep="")))
        eval(parse(text=paste("p_data_capx_q$",varname,"XlliqXsffr<-rlag(p_data_capx_q$",varname,",",as.character(lagc1_capx_q),")*rlag(p_data_capx_q$lliq_ind,",as.character(lagc1_capx_q),")*rlag(p_data_capx_q$sffr_wd,",as.character(lag0_capx_q),")",sep="")))
    } else {
        eval(parse(text=paste("p_data_capx_q$",varname,"Xsffr<-rlag(p_data_capx_q$",varname,",",as.character(lagc1_capx_q),")*rlag(p_data_capx_q$sffr,",as.character(lag0_capx_q),")",sep="")))
        eval(parse(text=paste("p_data_capx_q$",varname,"XlliqXsffr<-rlag(p_data_capx_q$",varname,",",as.character(lagc1_capx_q),")*rlag(p_data_capx_q$lliq_ind,",as.character(lagc1_capx_q),")*rlag(p_data_capx_q$sffr,",as.character(lag0_capx_q),")",sep="")))
    }
}
# Also, create the "selector_q" i.e. tobin_q variable to be instrumented, and interacted with "lliq_ind"
eval(parse(text=paste("p_data_capx_q$endg_var <-rlag(p_data_capx_q$",selector_q,",",as.character(h-h_q),")",sep="")))
eval(parse(text=paste("p_data_capx_q$endg_varXlliq <-rlag(p_data_capx_q$",selector_q,",",as.character(h-h_q),")*rlag(p_data_capx_q$lliq_ind,",as.character(h+1),")",sep="")))

### Non-cross-terms
# Add tover (fi_inst), lagged dependent variable, and other controls if requested
if (ALLin_ind){
    # Firm controls
    controls_nct_capx_q = c(paste("rlag(",fi_inst,", lagc1_capx_q)", sep=""), "rlag(lliq_ind, lagc1_capx_q)", paste("rlag(",fi_inst,", lagc1_capx_q):rlag(lliq_ind, lagc1_capx_q)", sep=""), firm_contrs, paste("rlag(lliq_ind, lagc1_capx_q):(", paste(firm_contrs, collapse=" + "), ")", sep=""), paste("rlag(",dep_var,", lagc1_capx_q)", sep=""), paste("rlag(lliq_ind, lagc1_capx_q):rlag(",dep_var,", lagc1_capx_q)", sep=""), "rlag(lev_rat, lagc1_capx_q)", "rlag(latq, lagc1_capx_q)", "rlag(cheat_rat, lagc1_capx_q)", "rlag(lliq_ind, lagc1_capx_q):rlag(lev_rat, lagc1_capx_q)", "rlag(lliq_ind, lagc1_capx_q):rlag(latq, lagc1_capx_q)", "rlag(lliq_ind, lagc1_capx_q):rlag(cheat_rat, lagc1_capx_q)")
} else {
    if (!fc_ind) {
        # No firm controls
        controls_nct_capx_q =c(paste("rlag(",fi_inst,", lagc1_capx_q)", sep=""), "rlag(lliq_ind, lagc1_capx_q)", paste("rlag(",fi_inst,", lagc1_capx_q):rlag(lliq_ind, lagc1_capx_q)", sep=""), paste("rlag(",dep_var,", lagc1_capx_q)", sep=""), paste("rlag(lliq_ind, lagc1_capx_q):rlag(",dep_var,", lagc1_capx_q)", sep=""))
    } else {
        # With firm controls
        controls_nct_capx_q =c(paste("rlag(",fi_inst,", lagc1_capx_q)", sep=""), "rlag(lliq_ind, lagc1_capx_q)", paste("rlag(",fi_inst,", lagc1_capx_q):rlag(lliq_ind, lagc1_capx_q)", sep=""), paste("rlag(",dep_var,", lagc1_capx_q)", sep=""), paste("rlag(lliq_ind, lagc1_capx_q):rlag(",dep_var,", lagc1_capx_q)", sep=""), "rlag(lev_rat, lagc1_capx_q)", "rlag(latq, lagc1_capx_q)", "rlag(cheat_rat, lagc1_capx_q)", "rlag(lliq_ind, lagc1_capx_q):rlag(lev_rat, lagc1_capx_q)", "rlag(lliq_ind, lagc1_capx_q):rlag(latq, lagc1_capx_q)", "rlag(lliq_ind, lagc1_capx_q):rlag(cheat_rat, lagc1_capx_q)")
    }
}

# Set up LHS variable and panel for running regression
p_data_capx_q["diff_dep_var"] <- p_data_capx_q[dep_var]
p_data_capx_q_used <- p_data_capx_q

### Control dummies
# Create correctly lagged (sic2_datacqtr X lliq) FE
p_data_capx_q_used$sic2_datacqtrXlliq_cur <- rlag(p_data_capx_q_used$sic2_datacqtrXlliq,lagc1_capx_q) # SIC2 X time X lliq FE
if (indt_ind) {
    dummies_capx_q = c("sic2_datacqtrXlliq_cur")
} else {
    if (ind_ind) {
        # Or if not, and only industry dummies requested
        dummies_capx_q = c("datacqtr", "sic1")
    } else {
        # Otherwise no industry dummies
        dummies_capx_q = c("datacqtr")
    }
}
# Create correctly lagged (firm X lliq) FE
p_data_capx_q_used$gvkeyXlliq_cur <- rlag(p_data_capx_q_used$gvkeyXlliq,lagc1_capx_q) # Firm X lliq FE
if (fe_ind) {
    dummies_capx_q = c(dummies_capx_q, "gvkeyXlliq_cur")
}

# Formula
if (cluftn_ind){
    formla_capx_q <- formula(paste("diff_dep_var ~ ", paste(controls_nct_capx_q, collapse=" + "), " + rlag(sffr",shock_name_ext,", lag0_capx_q)*(", paste(controls_ct_capx_q, collapse=" + "), ")| ", paste(dummies_capx_q, collapse=" + "), "|(endg_var + endg_varXlliq ~", paste(paste(fi_inst, "Xsffr",sep=""), collapse="+"),"+", paste(paste(fi_inst, "XlliqXsffr",sep=""), collapse="+"),") | sic3_datacqtr + gvkey" ))
} else{
    if (cluft_ind) {
        formla_capx_q <- formula(paste("diff_dep_var ~ ", paste(controls_nct_capx_q, collapse=" + "), " + rlag(sffr",shock_name_ext,", lag0_capx_q)*(", paste(controls_ct_capx_q, collapse=" + "), ")| ", paste(dummies_capx_q, collapse=" + "), "|(endg_var + endg_varXlliq ~", paste(paste(fi_inst, "Xsffr",sep=""), collapse="+"),"+", paste(paste(fi_inst, "XlliqXsffr",sep=""), collapse="+"),") | datacqtr + gvkey" ))
    } else {
        if (clu_ind){
            formla_capx_q <- formula(paste("diff_dep_var ~ ", paste(controls_nct_capx_q, collapse=" + "), "  | ", paste(dummies_capx_q, collapse=" + "), "|(endg_var + endg_varXlliq ~", paste(paste(fi_inst, "Xsffr",sep=""), collapse="+"),"+", paste(paste(fi_inst, "XlliqXsffr",sep=""), collapse="+"),") | datacqtr" ))
        } else{
            formla_capx_q <- formula(paste("diff_dep_var ~ ", paste(controls_nct_capx_q, collapse=" + "), "  | ", paste(dummies_capx_q, collapse=" + "), "|(endg_var + endg_varXlliq ~", paste(paste(fi_inst, "Xsffr",sep=""), collapse="+"),"+", paste(paste(fi_inst, "XlliqXsffr",sep=""), collapse="+"),")" ))
        }
    }
}

# Estimate
reg_capx_q <- felm(formla_capx_q, data=p_data_capx_q_used, na.action=na.exclude)
reg_sum_capx_q <- summary(reg_capx_q)
message(print(reg_sum_capx_q))

reg_coefs_capx_q <- coef(reg_sum_capx_q)
# Save coefs
betas_vec_q[ctr,1] <- reg_coefs_capx_q[, "Estimate"][paste("`endg_var(fit)`", sep="")]
ci_vec_q[ctr,] <- confint(reg_capx_q, level=ci_lev_set)[c(paste("`endg_var(fit)`", sep="")),]
# Save coefs for lliq_ind
yb_vec_q[ctr,1] <- reg_coefs_capx_q[, "Estimate"][paste("`endg_varXlliq(fit)`", sep="")]
yb_ci_vec_q[ctr,] <- confint(reg_capx_q, level=ci_lev_set)[c(paste("`endg_varXlliq(fit)`", sep="")),]

# Compute and collect correct ("true") confidence intervals of sums of parameters
if (ci_true_ind){
    crit_val = 1.96
    if (ci_lev_set==0.90){crit_val = 1.645}
    VCV_cur = reg_capx_q$clustervcv
    nK = dim(VCV_cur)[1]
    sel_vec <- numeric(nK)
    lev_ind <- match("`endg_var(fit)`", row.names(VCV_cur))
    yb_ind  <- match("`endg_varXlliq(fit)`", row.names(VCV_cur))
    sel_vec[c(lev_ind,yb_ind)] <- c(1,1)
    SE_sum <- sqrt(sel_vec %*% VCV_cur %*% sel_vec)
    # Save
    ylevb_vec_q[ctr,1] <- betas_vec_q[ctr] + yb_vec_q[ctr]
    ylevb_ci_vec_q[ctr,] <- ylevb_vec_q[ctr,1] + c(-1,1) * crit_val * SE_sum
}

} # End of estimations loop


# Create folder for figures
if (!dir.exists(paste("./saved_figs/",fig_name,sep=""))) {dir.create(paste("./saved_figs/",fig_name,sep=""))}

# Plot results for main coefs and coefs + lliq_ind
if (ci_true_ind){
    pdf(paste("./saved_figs/",fig_name,"/slevs_betas-yb_vec_q_citrue(",dep_var,")(",selector_q,")(i-",fi_inst,")",str_end,".pdf",sep=""), width=8, height=7)
    par(mar = c(5,5.3,2,2))
    cl<- rainbow(1, start=0.7)
    cl <- rgb(0.7, 0.7, 0.9)
    cl2 <- c("red")
    markers <- rep(c(4, 19), 1)
    if (ylim_man_ind){
        plot(lag_values, betas_vec_q[,1], las=1, type="n", ylim=ylim_man, xlab = "Quarters (h)", ylab=expression(paste(gamma["h"],", ", gamma["h"]+tilde(gamma)["h"])), cex.lab=font_lab_val, cex.axis=font_ax_val)
    } else {
        plot(lag_values, betas_vec_q[,1], type="n", ylim=c(-0.0002+min(0.00, min(ci_vec_q[,1], ylevb_ci_vec_q[,1])) , 0.0002+max(0,max(ci_vec_q[,2], ylevb_ci_vec_q[,2]))), xlab = "Quarters (h)", ylab=expression(paste(gamma["h"],", ", gamma["h"]+tilde(gamma)["h"])), cex.lab=font_lab_val, cex.axis=font_ax_val)
    }
    lines(lag_values,betas_vec_q[,1], type="o", col=cl[1], pch=markers[1], lwd=linew_main); lines(lag_values, ci_vec_q[,1], lty=2, col=cl[1], lwd=linew_ci); lines(lag_values,ci_vec_q[,2], lty=2, col=cl[1], lwd=linew_ci);
    lines(lag_values,ylevb_vec_q[,1], type="o", col=cl2[1], pch=markers[2], lwd=linew_main); lines(lag_values, ylevb_ci_vec_q[,1], lty=2, col=cl2[1], lwd=linew_ci); lines(lag_values, ylevb_ci_vec_q[,2], lty=2, col=cl2[1], lwd=linew_ci);
    abline(h=0.0); grid (NULL,NULL, lty = 6, col = "cornsilk2"); legend("bottomleft", lwd=rep(linew_main,2), lty=rep(2,1), col=c(cl[1],cl2[1]), pch=markers, cex=font_leg_val, legend=c(hliq_label,lliq_label))
    polygon(c(lag_values, rev(lag_values)), c(ylevb_ci_vec_q[,2], rev(ylevb_ci_vec_q[,1])), col = rgb(1,0,0, alpha=0.05), lty = 0)
    dev.off()
}

} # End of "dep_var" loop

# If ran Figure 4 for linv_rat_realh, then save the h=4 coefficients for back-of-the-envelope calculations in Section 5, and save the estimation results into a csv-file, to be then used by the structural model code when comparing the two
if ((fig_name %in% c("_fig4")) & (dep_var=="linv_rat_realh")){
    gamma_h_linv_Hliq = betas_vec_q[5]
    gamma_h_linv_Lliq = betas_vec_q[5] + yb_vec_q[5]

    # Save into csv-files
    if (TRUE){
        write.csv(betas_vec_q, paste("./saved_figs/slevsbetas_vec_q_vals(",dep_var,")(",selector_q,")",str_end,".csv",sep=""))
        write.csv(ci_vec_q, paste("./saved_figs/slevsbetas_vec_q_cis(",dep_var,")(",selector_q,")",str_end,".csv",sep=""))
        write.csv(ylevb_vec_q, paste("./saved_figs/slevsbetas_ylevb_vec_q_vals(",dep_var,")(",selector_q,")",str_end,".csv",sep=""))
        write.csv(ylevb_ci_vec_q, paste("./saved_figs/slevsbetas_ylevb_vec_q_cis(",dep_var,")(",selector_q,")",str_end,".csv",sep=""))
    }
}

# If ran Figure 13 and regressions are done, set sffr back to baseline
if (fig_name %in% c("_fig13")){
    p_data_capx_q$sffr <- p_data_capx_q$sffr_ff4 / sd(aggregate(sffr_ff4 ~ cqtr_num, data[(data$cqtr_num>=1990) & (data$cqtr_num<=2017),], function(x) mean(x))[,'sffr_ff4'])
}

} # fig_name loop end


######
# Construct Table 2
######
# Split across high/low liquidity
library(xtable)
sum_var_list <- c("atq_real", "age", "lev_rat", "cheat_rat", "mtover", "dylsaleq_real", "tobin_q", "gdp_beta", "MktRF_beta", "SMB_beta", "HML_beta", "RETXsd")
tabsave_rownames <- c("Total assets (2009 $MM)", "Age (years)", "Leverage (%)", "Liquidity ratio (%)", "Turnover (x100)", "Annual sales growth (%)", "Tobin's q", "GDP beta", "Market beta", "SMB beta", "HML beta", "Return volatility (%)")

dtl_mtover     <- sapply(p_data_capx_q[, sum_var_list], median, na.rm=TRUE)
dtl_H_mtover   <- sapply(p_data_capx_q[p_data_capx_q$Hmtover_ind==1, sum_var_list], median, na.rm=TRUE)
dtl_L_mtover   <- sapply(p_data_capx_q[p_data_capx_q$Hmtover_ind==0, sum_var_list], median, na.rm=TRUE)
# Combine into table
dtf_mtover     <- data.frame(dtl_mtover, dtl_H_mtover, dtl_L_mtover)
dtf_mtover[c("mtover"),] <- dtf_mtover[c("mtover"),]*avg_sdmtover*100
dtf_mtover[c("dylsaleq_real"),] <- dtf_mtover[c("dylsaleq_real"),]*100
colnames(dtf_mtover) <- c("All firms", "High turnover", "Low turnover")
rownames(dtf_mtover) <- tabsave_rownames
# Table for latex
xt_dtf_mtover <- xtable(dtf_mtover)
align(xt_dtf_mtover) <- c('l','c','c','c')
print(xt_dtf_mtover)


######
# Estimate regression (27)
######
formla_capx_q <- formula("ltobin_q ~ rlag(mtover, 1) + rlag(ltobin_q, 1) + rlag(sffr, 0) + rlag(mtover, 1):rlag(sffr, 0) | gvkey|0|0")
# Estimate
reg_capx_q <- felm(formla_capx_q, data=p_data_capx_q, na.action=na.exclude)
summary(reg_capx_q)
# Save coefficients
beta_0_ntd  = as.double(reg_capx_q$coefficients["rlag(sffr, 0)",])
gamma_0_ntd = as.double(reg_capx_q$coefficients["rlag(mtover, 1):rlag(sffr, 0)",])
# Impose coefficients manually
# beta_0_ntd  = -0.3858361
# gamma_0_ntd = -0.0974441


######
# Back-of-the-envelope calculations in Section 5
######
# Exposure of Fed funds rate change to sffr
# Back out aggregate time series from panel
library(data.table)
agg_sffr_ff4 = data.table(p_data_capx_q)[, .(sffr_ff4 = mean(sffr_ff4)), by = .(cqtr_num)]
agg_dffr     = data.table(p_data_capx_q)[, .(dffr = mean(dffr)), by = .(cqtr_num)]
# Estimate
reg_agg <- lm(agg_dffr$dffr ~ agg_sffr_ff4$sffr_ff4)
# Safe coefficient
dffr_dsffr  = as.double(reg_agg$coefficients["agg_sffr_ff4$sffr_ff4"])
# dffr_dsffr  = 2.97818
sffr_sd     = sd(aggregate(sffr_ff4 ~ cqtr_num, data[(data$cqtr_num>=1990) & (data$cqtr_num<=2017),], function(x) mean(x))[,'sffr_ff4'])
# sffr_sd     = 0.096621979654003
# Share of public firms in aggregate investment (Asker et al., 2012)
pf_share    = 1-0.545

# IV coefficient values (log(inv) on q) imposed already above in estimation
# gamma_h_linv_Hliq = -0.34857
# gamma_h_linv_Lliq = gamma_h_linv_Hliq + 1.50679

# No time dummies' effect on q
p_data_capx_q["gammaXmtover_ntd"] = ( beta_0_ntd  + gamma_0_ntd * p_data_capx_q$mtover)  / (100*sffr_sd) # Per basis point monetary shock eps^m
# And for interaction with gamma_h from investment rate IV
p_data_capx_q["lghgammaXmtover_ntd"] = ( (1-p_data_capx_q$lliq_ind)*gamma_h_linv_Hliq + p_data_capx_q$lliq_ind*gamma_h_linv_Lliq ) * p_data_capx_q$gammaXmtover_ntd # "Per basis point monetary shock"

# Histogram of stock price semi-elasticities (Figure 5)
if (!dir.exists("./saved_figs/_fig6")) {dir.create("./saved_figs/_fig6")}
pdf(paste("./saved_figs/_fig6/hist_lghgammaXmtover_ntd.pdf",sep=""), width=8, height=7)
hist(p_data_capx_q[,"lghgammaXmtover_ntd"], xlab=TeX('$\\left(\\gamma_{4} + \\tilde{\\gamma}_{4} \\mathbf{I}_{t-1}^{i} \\right) \\left(\\beta_{2} + \\beta_{3} \\tau_{t-1}^{i} \\right)$'), breaks=60, freq=TRUE, main=NULL, xlim=c(-0.14,0.05), axes=FALSE, col="white")
axis(1,pos=0)
axis(2,pos=-0.145)
dev.off()


# Compute weighted avgs
# Precreate a "weight capxq" variable to equal 0 when NA
p_data_capx_q["wcapxq"] <- p_data_capx_q["capxq"]
p_data_capx_q[is.na(p_data_capx_q$capxq), "wcapxq"] <- 0.0
# Weighted aggregate by quarter
aggpf_ilinv_resp_ntd = mean(data.table(p_data_capx_q)[, list(wmt = weighted.mean(lghgammaXmtover_ntd, wcapxq, na.rm=TRUE)), by=cqtr_num]$wmt, na.rm=TRUE)
# Report correctly scaled aggregates:
100* sffr_sd * aggpf_ilinv_resp_ntd  # Aggregated public firm investment response (in percentages) to 1 sd sffr (9.66bp) increase (reported first in Section 5)
100* 0.25*aggpf_ilinv_resp_ntd       # Aggregated public firm investment response (in percentages) to 25bp sffr increase (reported second in Section 5)
100* 0.25*aggpf_ilinv_resp_ntd / dffr_dsffr # Aggregated public firm investment response (in percentages) to 25bp sffr increase (reported third in Section 5)
# Taking into account public firms' share in the aggregate
100* 0.25*aggpf_ilinv_resp_ntd * pf_share / dffr_dsffr # Aggregated firm investment response (in percentages) to 25bp sffr increase (reported fourth in Section 5)


# Compute the internal calibration target moments (Table 1)
if (TRUE){
    # 1. Median cash-to-assets (0.086)
    round(median(p_data_capx_q$cheat_rat, na.rm=TRUE)/100, 3)
    # 2. Median I/K for low cash-to-assets firms (0.039)
    round(median(p_data_capx_q[(rlag(p_data_capx_q$lliq_ind,1)==1), "inv_rat_realh"], na.rm=TRUE)/100, 3)
    # 3. Median I/K for high cash-to-assets firms (0.056)
    round(median(p_data_capx_q[(rlag(p_data_capx_q$lliq_ind,1)==0), "inv_rat_realh"], na.rm=TRUE)/100, 3)
    # 4. Frequency of equity issuance (0.080)
    round(sum((p_data_capx_q$neqat_rat_realh > 5/4), na.rm=TRUE) / sum(!(is.na(p_data_capx_q[,"neqat_rat_realh"]))), 3)
    # 5. Mean issuance conditional on issuing (0.157):
    round(mean(p_data_capx_q[(p_data_capx_q$neqat_rat_realh > 5/4),"neqat_rat_realh"], na.rm=TRUE)/100, 3)
}
